#include <math.h>

#include <iostream>

#include "Interpolation3D.h"
#include "test.h"

typedef double (*Func)(double, double, double);

template <class InternalType = double, class ExternalType = float>
void test_subcase(Interpolation3DBase<InternalType, ExternalType> &b,
                  ExternalType xa[], ExternalType ya[], ExternalType za[],
                  ExternalType ***va, int m, int n, int l, Func f) {
  b.construct(xa, ya, za, va);
  for (int i = 0; i < m; i++) {
    for (int j = 0; j < n; j++) {
      for (int k = 0; k < l; k++) {
        ExternalType x = sqrt(i);
        ExternalType y = sqrt(j);
        ExternalType z = sqrt(k);
        // 1e-6 because float has 6 exact numbers
        ExternalType error =
            (f(x, y, z) == 0)
                ? fabs(b.evaluate(x, y, z) - f(x, y, z))
                : fabs(b.evaluate(x, y, z) - f(x, y, z)) / f(x, y, z);
        CHECK(error < 1e-6);
        // printf("x = %.1f, y = %.1f, f = %.2f\n", x, y, b.evaluate(x, y));
      }
    }
  }
}

template <class InternalType = double, class ExternalType = float>
void test_case(int m, int n, int l, Func f, int deg) {
  ExternalType *xa = new ExternalType[m];
  ExternalType *ya = new ExternalType[n];
  ExternalType *za = new ExternalType[l];
  ExternalType ***va = new ExternalType **[m];
  for (int i = 0; i < m; i++) xa[i] = sqrt(i);
  for (int j = 0; j < n; j++) ya[j] = sqrt(j);
  for (int k = 0; k < l; k++) za[k] = sqrt(k);
  for (int i = 0; i < m; i++) {
    va[i] = new ExternalType *[n];
    for (int j = 0; j < n; j++) {
      va[i][j] = new ExternalType[l];
      for (int k = 0; k < l; k++) va[i][j][k] = f(xa[i], ya[j], za[k]);
    }
  }
  SUBCASE("Efficient Linear Interpolation") {
    // std::cout << "Efficient Linear Interpolation" << std::endl;
    LinInt3D<InternalType, ExternalType> p =
        LinInt3D<InternalType, ExternalType>(m, n, l);
    test_subcase(p, xa, ya, za, va, m, n, l, f);
  }
  SUBCASE("Quadratic Interpolation") {
    // std::cout << "Quadratic Interpolation" << std::endl;
    QuaInt3D<InternalType, ExternalType> h =
        QuaInt3D<InternalType, ExternalType>(m, n, l);
    test_subcase(h, xa, ya, za, va, m, n, l, f);
  }
  SUBCASE("Spline Interpolation") {
    // std::cout << "Spline Interpolation" << std::endl;
    SplInt3D<InternalType, ExternalType> s =
        SplInt3D<InternalType, ExternalType>(m, n, l);
    test_subcase(s, xa, ya, za, va, m, n, l, f);
  }
  for (int i = 0; i < m; i++)
    for (int j = 0; j < n; j++) delete[] va[i][j];
  for (int i = 0; i < m; i++) delete[] va[i];
  delete[] xa, ya, za;
}

TEST_CASE("Test: Linear Polynomial") {
  int m = 10, n = 10, l = 10;
  test_case<double, float>(m, n, l, plinear, 1);
}

// TEST_CASE("Test: Bi-Quadratic Polynomial") {
//   int m = 10, n = 10;
//   test_case<double, float>(m, n, pquadratic, 2);
// }

// TEST_CASE("Test: Bi-Cubic Polynomial") {
//   int m = 10, n = 10;
//   test_case<double, float>(m, n, pcubic, 3);
// }

TEST_CASE("Test: Trigonometric Function") {
  int m = 10, n = 10, l = 10;
  test_case<double, float>(m, n, l, fsin, 3);
}
